Demonstration of the IBCAO class

Plot map with some data on it

Load necessary Python modules



In [1]:

    
# This is designed to be run from the doc/ directory
import sys
sys.path.append ('..')



In [2]:

    
# load matplotlib and cartopy
%matplotlib notebook
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import cartopy.crs as ccrs
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER

# load the IBCAO class
from ibcao import *









    



/usr/lib/python3.6/site-packages/matplotlib/ticker.py:1685: UserWarning: Steps argument should be a sequence of numbers
increasing from 1 to 10, inclusive. Behavior with
values outside this range is undefined, and will
raise a ValueError in future versions of mpl.
  warnings.warn('Steps argument should be a sequence of numbers\n'

Set up the basic IBCAO figure



In [3]:

    
i = IBCAO ()
f = i.template ()
f.set_size_inches (10,10)

## We will get a warning about an open netcdf file, this is caused by the map data being memory-mapped
## so that only the necessary data is loaded. See https://github.com/gauteh/ibcao_py/issues/1 for more
## information.









    



ibcao read, shape: (11617, 11617)

Let us add some longitude and latitude markers to the map



In [4]:

    
f = i.template ()
f.set_size_inches (10,10)

# lets put some text along the parallels
lat = np.arange (65, 90, 5)
lon = np.repeat (0, len(lat))

# regular lat, lon projection
g = ccrs.Geodetic ()

for lon, lat in zip (lon, lat):
    plt.text (lon, lat, LATITUDE_FORMATTER.format_data(lat), transform = g)

# and some text along the meridians
lon = [-45, 45, 135, -135]
lat = np.repeat (70, len(lon))

for lon, lat in zip (lon, lat):
    plt.text (lon, lat, LONGITUDE_FORMATTER.format_data(lon), transform = g)

# also; the north pole
plt.text (0, 90, "NP", transform = g)









    














    











    Out[4]:





<matplotlib.text.Text at 0x7f56f1d84630>

Plotting data on the map

If we have a data set, given in latitude and longitude degrees, we can plot this directly using the Geodetic transform.



In [5]:

    
data = np.array( [[75, 35, 10],
                  [76, 34, 9],
                  [80, 17, 11],
                  [85, 15, 20]]) 
f = i.template ()
f.set_size_inches (10, 10)

## Note that longitude = x, latitude = y
plt.scatter (data[:,1], data[:,0], s = 10 * data[:,2], c = np.random.rand(4), transform = i.g)
plt.plot (data[:,1], data[:,0], transform = i.g, linewidth = 2)

## Show a smaller section of the map (coordinates given in degrees)
f.gca().set_extent ((0, 45, 70, 85))

Get depth along intersect

Get and plot depth profile along the great circle between two points, say the Fram strait between Svalbard and Greenland.



In [6]:

    
f = i.template ()
f.set_size_inches (10,10)

## Zoom in on the Fram Strait
f.gca().set_extent ((-18, 16, 75, 81))

## Start and end points
start = [10, 78]
end   = [-18, 76]
plt.plot (start[0], start[1], 'o', markersize = 20, transform = i.g)
plt.plot (end[0], end[1], 'o', markersize = 20, transform = i.g)

## Sample a great circle between the two points using the pyproj.Geod set up by the IBCAO class
gc = np.array(i.geod.npts (start[0], start[1], end[0], end[1], 100))
plt.plot (gc[:,0], gc[:,1], transform = i.g, linewidth = 2)

## Interpolate the depth along the great circle
gc_xy = i.ups.transform_points (i.g, gc[:,0], gc[:,1]) # convert to UPS coordinates
depth = i.map_depth (gc_xy[:,0], gc_xy[:,1])

## Plot depth profile along latitude
df = plt.figure (figsize = (10,5))
plt.plot (gc[:,0], depth)
plt.xlabel ('Longitude [deg]')
plt.ylabel ('Depth [m]')









    














    











    














    











    Out[6]:





<matplotlib.text.Text at 0x7f56b152fef0>